rm(list=ls())
setwd("/Users/afanlo/Documents/Projects/Kroenig Rep")
library(dplyr)
library(MASS)
library(xtable)
library(stargazer)
library(lmtest)
library(sandwich)


### Creating a function to compute clustered standard errors
vcovCluster <- function(
    model,
    cluster
)
{
  require(sandwich)
  require(lmtest)
  
  cluster = as.factor(cluster)
  
  if (!is.null(model$na.action)){
    omit.rows = model$na.action
    cluster = cluster[-omit.rows]
  }
  if (nrow(model.matrix(model)) != length(cluster)){
    stop("something's not working: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap or randomization inference instead).")
  }
  
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- na.omit(apply(estfun(model), 2, function(x) tapply(x, cluster, sum)))
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##
## New Superiority Measure with Original ICB Outcome Codings
##

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))


##Changing back to ICB coding for Korean War
df$victory<-ifelse(df$crisno==133 & df$ccode==365,  1, df$victory)

## Changing back to ICB Coding of 1969 Sino-Soviet Border War
df$victory<-ifelse(df$crisno==231, 0, df$victory)


## NEW SUPERIORITY MEASURE
df$sup_new<-ifelse(df$nucratio>3 & df$strikea==1, 1, 0)



## Models for Original ICB Outcome Codings without Stakes
mod1<-glm(victory~sup_new + proximity + polity21 + 
            capshare + strikea +  tpop_1 + viol + cris_ave, data=df, family=binomial(link="logit"))
mod1test<-coeftest(mod1, vcov=vcovCluster(model=mod1, cluster=factor(df$crisisdyad)))

## Models for Original ICB Outcome Codings with Stakes
mod2<-glm(victory~sup_new + stakes + proximity + polity21 + 
            capshare + strikea +  tpop_1 + viol + cris_ave, data=df, family=binomial(link="logit"))
mod2test<-coeftest(mod2, vcov=vcovCluster(model=mod2, cluster=factor(df$crisisdyad)))

stargazer(mod1test, mod2test)

##
## New Superiority Measure with Kroenig's Outcome Codings
##

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))


## NEW SUPERIORITY MEASURE
df$sup_new<-ifelse(df$nucratio>3 & df$strikea==1, 1, 0)

## Models for Kroenig's Outcome Codings without Stakes
mod3<-glm(victory~sup_new + proximity + polity21 + 
            capshare + strikea +  tpop_1 + viol + cris_ave, data=df, family=binomial(link="logit"))
mod3test<-coeftest(mod3, vcov=vcovCluster(model=mod3, cluster=factor(df$crisisdyad)))

## Models for Kroenig's Outcome Codings with Stakes
mod4<-glm(victory~sup_new + stakes + proximity + polity21 + 
            capshare + strikea +  tpop_1 + viol + cris_ave, data=df, family=binomial(link="logit"))
mod4test<-coeftest(mod4, vcov=vcovCluster(model=mod4, cluster=factor(df$crisisdyad)))


## Table 6
stargazer(mod1test, mod2test, mod3test, mod4test)

##
## RANDOMIZATION INFERENCE (Added to Table 6, Columns 1 and 3)
##

rm(list=ls())
set.seed(02303)
library(dplyr)
library(MASS)
library(xtable)
library(stargazer)
library(lmtest)
library(sandwich)
###
##Running randomization inference with the 3x superiority measure and Original ICB Outcome Codings. 
###


setwd("/Users/afanlo/Documents/Projects/Kroenig Rep")

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))

##Changing back to ICB coding for Korean War
df$victory<-ifelse(df$crisno==133 & df$ccode==365,  1, df$victory)

## Changing back to ICB Coding of 1969 Sino-Soviet Border War
df$victory<-ifelse(df$crisno==231, 0, df$victory)

## NEW SUPERIORITY MEASURE
df$sup_new<-ifelse(df$nucratio>3 & df$strikea==1, 1, 0)

bin_df3x<-glm(data=df, victory~sup_new+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
sims<-10000

coef.bindf3x.ss<-matrix(NA, sims, 9)
coefs<-matrix(NA, sims, 9)
set.seed(1e-08)

reps<-sum(df$sup_new) ###NEW CODE
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has 3x arsenal of side b
    sidea<-c(rep(0,29-reps), rep(1, reps)) ###CHANGED
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 29, 2)
    c[,ran]<-sample(sidea, 29, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:58,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$sup_new_cut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(sup_new_cut))
    
    model<-glm(data=df.new, victory~sup_new_cut+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.bindf3x.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.bindf3x.ss<-cbind(coef.bindf3x.ss, fail)
coef.bindf3x.ss<-subset(coef.bindf3x.ss, coef.bindf3x.ss[,10]==0)

coef.bindf3x1.ss<-matrix(NA, 1, 9)
coef.bindf3x1.ss[1,]<-summary(bin_df3x)$coefficients[,3]

p.coef3x.1<-mean(sapply(coef.bindf3x.ss[,1], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,1])))
p.coef3x.2<-mean(sapply(coef.bindf3x.ss[,2], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,2])))
p.coef3x.3<-mean(sapply(coef.bindf3x.ss[,3], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,3])))
p.coef3x.4<-mean(sapply(coef.bindf3x.ss[,4], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,4])))
p.coef3x.5<-mean(sapply(coef.bindf3x.ss[,5], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,5])))
p.coef3x.6<-mean(sapply(coef.bindf3x.ss[,6], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,6])))
p.coef3x.7<-mean(sapply(coef.bindf3x.ss[,7], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,7])))
p.coef3x.8<-mean(sapply(coef.bindf3x.ss[,8], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,8])))
p.coef3x.9<-mean(sapply(coef.bindf3x.ss[,9], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,9])))

pvals_rand<-matrix(NA, 1, 9)
pvals_rand[1,]<-c(p.coef3x.1, p.coef3x.2, p.coef3x.3, p.coef3x.4, p.coef3x.5,
                  p.coef3x.6, p.coef3x.7, p.coef3x.8,  p.coef3x.9)

## Extract p-value for new superiority measure for original ICB codings and add to Table 6, Column 1
pvals_rand[,2]
###
#Running randomization inference with the  3x superiority measure and Kroenig Codings. 
###
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))


## NEW SUPERIORITY MEASURE
df$sup_new<-ifelse(df$nucratio>3 & df$strikea==1, 1, 0)

bin_df3x<-glm(data=df, victory~sup_new+proximity +  polity21 + 
                capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
sims<-10000

coef.bindf3x.ss<-matrix(NA, sims, 9)
coefs<-matrix(NA, sims, 9)
set.seed(1e-08)

reps<-sum(df$sup_new) ###NEW CODE
fail <- matrix(NA, 10000, 1)


for (i in 1:sims){
  tryCatch( expr={
    #Creating a random vector of zeros and ones to randomize whether side a has 3x arsenal of side b
    sidea<-c(rep(0,29-reps), rep(1, reps)) ###CHANGED
    #creating a matrix to store the randomization of side a- and then attaching zeros next to it so that I ensure the other country in the dyad doesn't also get that they would sustain less damage 
    ran<-as.numeric(sample(1:2, 1))
    other<-as.numeric(ifelse(ran==2, 1, 2))
    c<-matrix(NA, 29, 2)
    c[,ran]<-sample(sidea, 29, replace=F)
    c[,other]<-0
    c<-t(c)
    ass<-as.vector(c)
    #ordering the data frame by crisis and by dyad
    df.new<-df[1:58,]
    df.new<-df.new[order(df.new$crisisdyad),]
    
    #applying the assignment of the ordered side a variable
    df.new$sup_new_cut<-ass
    #making the ordered side b variable to be the opposite of the ordered side a variable
    df.new<-df.new %>% group_by(crisisdyad) %>% mutate(sum_ord= sum(sup_new_cut))
    
    model<-glm(data=df.new, victory~sup_new_cut+proximity +  polity21 + 
                 capshare + strikea +  tpop_1 + viol + cris_ave, family=binomial(link="logit"))
    coefs[i,]<-summary(model)$coefficients[,1]
    coef.bindf3x.ss[i,]<-summary(model)$coefficients[,3]
    df.new<-NA
    fail[i]<-0},
    
    warning = function(w) {
      print(w)
    }
  )}

fail[,1]<-ifelse(is.na(fail[,1])==T, 1, fail[,1])
coef.bindf3x.ss<-cbind(coef.bindf3x.ss, fail)
coef.bindf3x.ss<-subset(coef.bindf3x.ss, coef.bindf3x.ss[,10]==0)

coef.bindf3x1.ss<-matrix(NA, 1, 9)
coef.bindf3x1.ss[1,]<-summary(bin_df3x)$coefficients[,3]

p.coef3x.1<-mean(sapply(coef.bindf3x.ss[,1], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,1])))
p.coef3x.2<-mean(sapply(coef.bindf3x.ss[,2], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,2])))
p.coef3x.3<-mean(sapply(coef.bindf3x.ss[,3], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,3])))
p.coef3x.4<-mean(sapply(coef.bindf3x.ss[,4], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,4])))
p.coef3x.5<-mean(sapply(coef.bindf3x.ss[,5], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,5])))
p.coef3x.6<-mean(sapply(coef.bindf3x.ss[,6], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,6])))
p.coef3x.7<-mean(sapply(coef.bindf3x.ss[,7], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,7])))
p.coef3x.8<-mean(sapply(coef.bindf3x.ss[,8], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,8])))
p.coef3x.9<-mean(sapply(coef.bindf3x.ss[,9], FUN=function(x) x >= abs(coef.bindf3x1.ss[1,9])))

pvals_rand<-matrix(NA, 1, 9)
pvals_rand[1,]<-c(p.coef3x.1, p.coef3x.2, p.coef3x.3, p.coef3x.4, p.coef3x.5,
                  p.coef3x.6, p.coef3x.7, p.coef3x.8,  p.coef3x.9)


## Extract p-value for new superiority measure for Kroenig recodings and add to Table 6, Column 3
pvals_rand[,2]

